Domain regime in two-dimensional disordered vortex matter 
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A detailed numerical study of the real space configuration of vortices in disordered superconductors 
using 2D London-Langevin model is presented. The magnetic field B is varied between and Bc2 
for various pinning strengths A. For weak pinning, an inhomogeneous disordered vortex matter is 
observed, in which the topologically ordered vortex lattice survives in large domains. The majority of 
the dislocations in this state are confined to the grain boundaries/domain walls. Such quasi-ordered 
configurations are observed in the intermediate fields, and we refer it as the domain regime (DR). The 
DR is distinct from the low-field and the high-fields amorphous regimes which are characterized by a 
homogeneous distribution of defects over the entire system. Analysis of the real space configuration 
suggests domain wall roughening as a possible mechanism for the crossover from the DR to the 
high-field amorphous regime. The DR also shows a sharp crossover to the high temperature vortex 
hquid phase. The domain size distribution and the roughness exponent of the lattice in the DR are 
also calculated. The results are compared with some of the recent Bitter decoration experiments. 

PACS numbers: 74.25. Qt, 74.25. Sv 



INTRODUCTION 

The vortex state in type-II superconductors is a 
paradigm for studying the effect of quenched disorder in 
condensed matter. Over the last decade, much of the ef- 
fort has been spent on characterizing the various phases 
of the vortex state as a function of the magnetic field B 
and the temperature T. For 3Z? vortex systenyll, three 
phases have been identified unambiguously 0, Ull] '■ the 
Bragg glass (BG) with quasi-long range order, the amor- 
phous vortex glass (VG), and the vortex liquid (VL). The 
VG and the VL are distinguished by their superconduct- 
ing and ohmic responses, respectively. Experiments in 
high-Tc superconductors suggest that the BG phase ap- 
pears in the low-B and low-T region whereas the VG 
phase occupies the high-S and low-T region of the B- 
T phase diagram. The VL phase appears close to the 
upper critical field Bc2{T). 

The first detailed calculation of the real space 
structure of the vortex lattice in the presence of 
quenched impurities (sll was carried out by Larkin and 
Ovchinnikov[6| (LO). The LO theory assumes that for 
the weak pinning, the vortex lattice is coherently pinned 
within a volume Vc- Beyond Vc, the effect of impurities 
dominates and the long range positional order is lost. 
It was further proposed that when the vortex displace- 
ment becomes of the order of lattice constant, topologi- 
cal defects (dislocations in D — 2 and dislocation loops 
in _D = 3) are generated(7|. Later calculations by Gia- 
marchi and Le Doussal |3|,|3| (GLD) showed that the LO 
calculation overestimates the effect of impurities at large 
distances and the vortex displacement grows only loga- 
rithmically (for 31? vortex system). The positional cor- 
relation decays as C{r) ^ -pj, where 77 is a non-universal 
exponent'^. The quasi-long range G{r) leads to a topo- 
logically ordered phase (the Bragg glass), which is sta- 



ble with respect to the formation of defects[lO|. On in- 
creasing B (or the pinning strength), the LO and GLD 
theories predict proliferation of topological defects, thus 
forming the VG at high fields The transition be- 
tween the BG and the high field VG is predicted to be 
of first-order [H E H EB. The BG also undergoes a 
melting transition to VL on increasing T. 

In D = 2, the BG phase is unstable to the formation 
of dislocations and the positional quasi-long range order 
is destroyed 0, 0|. However, for weak pinning and at 
low temperatures, the unbound dislocations appear only 
at large length scale ^ Ra, where Ra is the "Random 
Manifold" length scale and is the distance at which po- 
sitional correlation begins to decayj^. On length scales 
shorter than (^jj, the topologically ordered lattice forms 
a quasi-Bragg glass (qBG). Such a qBG shows an ex- 
ponentially sharp crossover to the high temperature VL 
phase 01, reminiscent of the "melting" transition of the 
pure system. A similar exponential crossover was pro- 
posed between the qBG and the VG phase as a function 
of i?, or pinning strength. 

The real space structure of the vortex system has been 
studied used neutron diffraction Eol EH 12^ and Bitter 
decoration of vorticesjH HI El l2g|. The latter tech- 
nique allows direct visualization of the large scale struc- 
ture of the configuration and hence enables one to ana- 
lyze the role of topological defects on the decay of trans- 
lational order. Recent decoration experiments 
of NbSe2 have raised some important issues concerning 
the nature of the disordered phase. Previous transport 
measurements on the same samples of NbSe2 sug- 
gested an order-disorder transition on increasing B (or 
T). Fasano et al. showed that the spatial configura- 
tion of vortices does not show any significant difference 
between the ordered and the disordered vortex phases 
identified in Ref . [2^ . More importantly, both phases 
were found to be polycrystalline with dislocations form- 
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ing grain boundaries. Within each grain, the lattice 
shows significant bond orientational order. This is in 
contrast to the naive theoretical picture of the ordered 
phase which expects a dislocation-free configuration, and 
the disordered phase in which the distribution of the dis- 
locations is expected to be homogeneous. 

In this paper, we analyze in detail the real space con- 
figuration of the disordered phase using numerical sim- 
ulation of a 2D vortex system at T = 0. The magnetic 
field B is varied over a wide range for various values of 
the pinning strength A. The real space configuration 
shows that for the intermediate field range, the system 
shows inhomogeneous disordering. The majority of the 
dislocations are confined to the grain boundaries which 
forms the domain wall between regions of ordered lat- 
tice. The domain size and its distribution is dependent 
on B and A. We refer the intermediate fields in which 
the vortex state is quasi-ordered as the domain regime 
(DR) . The DR is distinct from the amorphous regime at 
low-fields and high-fields, where the defects appear at a 
length scale ~ Og (lattice constant) and its distribution 
is homogeneous over the entire system. Analysis of the 
real space images suggests domain wall roughening as 
a possible mechanism for the crossover between the DR 
and the high-fields amorphous regime. We also obtained 
the roughening exponent C, of the vortex lattice in the do- 
main regime. Finite temperature simulation shows that 
the domain regime undergoes a sharp crossover to the 
high-T liquid phase, which is reminiscent of the thermal 
melting in the pure vortex system. 

The paper is organized as follows: in section II, we 
discuss the simulation approach in detail. The results 
and the analysis of the real space configuration are pre- 
sented in section III, followed by conclusions in section 
IV. 

SIMULATION METHOD 

We consider a 2D cross-section perpendicular to the 
magnetic field B — Bi of a bulk type-II superconductor 
in the mixed state. Within London's approximation, the 
vortex can be considered as a point particle with the dy- 
namics governed by an overdamped equation of motion 

dr 

j^i k 

+Fe,t+CW. (1) 

Here, rj is the flux-flow viscosity. On the left hand 
side, the first term represents the inter-vortex interac- 
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tion U'"{r) = j^Ko{f/\), where ATo is the zeroth- 
order Bessel function, and f = (r^ -f 2^^)^/^. 00 is 
the flux quantum, and the A and ^ are the penetra- 
tion depth and the coherence length of the supercon- 
ductor, respectively. This form of the inter-vortex in- 



teraction includes the finite core size of the vortex j2q . 
The second term represents vortex pinning by parabolic 
potential wells, where [/^(r) = U[){^ — 1) for r < r„, 

and otherwise. The pinning centers are randomly lo- 
cated at positions in the simulation box. The third 
term Yext = ^ '/'o^ is the Lorentz force experienced 
by the vortex due to the transport current density J. 
The thermal noise is represented by C, with {Ci.p{t)) — Q, 
and {Ci,p{t)Cj,p' {t')) = 2kBTri6ijSpp,S{t - t'), where T is 
the temperature, ks is the Boltzmann constant, and 
p,p' — x,y. The length is in units of A(i? = 0, T = 0) = Ao, 
and the temperature T is in units of Ao/o/^b, where 
fo = g^2xi ■ The current density J and the velocity v of 
the vortices are in units of c/o/0o and fo/rj, respectively. 
Also, the prefactor for the pinning potential Uq is scaled 
by /o/Ao. 

We use the reduced magnetic field b = B/Bc2, where 
the upper critical field Bc2 = and £,o = £,{B = 

0,T = 0). The b is calculated from the lattice constant 
SI — (i^)2(--^)2. The Ginzburg-Landau parameter 

ii = J is an input to the simulation. The magnetic field 
dependence of the length scales A and ^ follows the rela- 
tion A(6) = /(fc)Ao and ^(6) = f{b)^Q, respectively. The 
renormalization factor f{b) = (1 — 6^)~2. This form 
of the renormalization factor is similar to the temper- 
ature dependence of ^ and A in the Ginzburg-Landau 
theoryHa with T/Tc replaced by [B / Bc2f . Similar 
form of the renormalization factor for A have been used 
in Ref. '30*1 . The parameters used in the simulation are 
n — IQ and Aq = 1000 A, which are close to the values for 
the low- Tc superconductors, particularly NbSe2. Peri- 
odic boundary conditions are imposed in both directions, 
and the minimum image convention is followed. The 
magnetic field b is varied by changing the size of the sim- 
ulation box, keeping the number of vortices A^„ — 4096 
fixed. Simulations were also performed using Ny be- 
tween 800-1200, and for some parameters A^„ — 6400 
was used to check for the finite size effects. The Uq is 
distributed randomly between AitO.Ol, where A = {Uq). 
The range of the pinning potential rp — ^q. In this pa- 
per, we present results for pin density Up = 2.315/Ao. 
For T = 0, Eq.(l) is time integrated by the predictor- 
corrector scheme, and the finite ternperature simulation 
is performed using Heun's method |3 If. The simulation 
at high vortex densities requires long computational time 
and parallel algorithms were employed to reduce the run 
time. Details of the implementation of the parallel algo- 
rithms can be found in Ref . . 

The real space configuration is characterized by the 
topological defect density n^/Ao (number of defects per 
unit area of the simulation box). Below, we also use 
the defect fraction /d, which is defined as the number of 
defects per vortex. The defects are defined as vortices 
with coordination number other than 6 and are identified 
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by Delaunay triangulation of the real space position of 
the vortices. In 2D systems, the vortices with coordina- 
tion number 5 and 7 are dischnations. A 5-dischnation 
and a 7-dischnation separated by a distance oq forms a 
bound pair which is an edge dislocation. Over most of 
the field range, the fraction of free dischnations is neg- 
ligibly small and the majority of the defects are edge 
dislocations. Hence, the defect density Ud is approxi- 
mately twice the dislocation density in the system. The 
hexatic order in the system is quantified by the six-fold 
orientational order parameter '^q — e^^'^l, where 

6ij is the angle between the nearest neighbor vortices 
relative to a reference axis. 

The simulation is performed by two different meth- 
ods. In the first method, we start with a perfect vor- 
tex lattice and the driving current /(oc J) is reduced 
to from a value much greater than the depinning cur- 
rent Ic- This is referred as the current annealing (CA) 
method. In the second method, the conventional ther- 
mal annealing (TA) is applied wherein the temperature 
T is reduced to in small steps from the high temper- 
ature liquid phase (also known as simulated annealing). 
Experimentally, the TA is equivalent to the field cool- 
ing procedure. We have shown previously|3^ that the 
configuration obtained by CA is stable to small pertur- 
bations compared to the configuration obtained by TA. 
This is also supported by experiments j^li which show 
that the field cooled state is unstable to small driving 
force I <^ Ic and a stable configuration is obtained when 
the system is brought to rest after driven with I ^ Ic- 
The two methods, CA and TA, are compared in section 
III. B. 



RESULTS AND DISCUSSIONS 
Zero temperature simulation 

In this section, we analyze the zero-temperature con- 
figurations obtained by the current annealing method. 
The system is slowly brought to rest across the depin- 
ning current for each value of the magnetic field b. In 
the absence of thermal fluctuations, the vortex config- 
uration is determined by the balance between the long 
range elastic force and the pinning force. We first show 
the real space images of the configuration as the mag- 
netic field b is increased. 

Real space configuration 

Fig.l shows the Delaunay triangulation in a region of 
the simulation box for various values of the magnetic 
field. The pinning strength A = 0.02, and iV„ = 4096, 
except for b = 0.1 for which N^^ = 900. At small fields 
b < 0.1, the defect distribution is homogeneous over the 




FIG. 1: (Color Online) The real space configuration of vor- 
tices in a region of the simulation box. The values of b are 
(a) 0.1, (b) 0.4, (c) 0.5, (d) 0.6, (e) 0.8, and (f) 0.9. The 
black (red) and gray (blue) dots denote vortices with 7 and 
5 neighbors, respectively j^a]- For b = 0.1, = 900, and 
for the rest of the images A'^,, — 4096. The pinning strength 
A = 0.02. 

entire system and the configuration is amorphous. The 
defect fraction (number of defects per vortex) fd > 0.35 
at these low fields. 

With increasing 6, small regions of ordered lattice start 
appearing. This can be seen for 6 ~ 0.1 in which ordered 
lattice is formed in regions less than 3-4ao wide. For 
b > 0.2, the defect distribution becomes inhomogeneous. 
The dislocations come closer to form a network of grain 
boundaries across the system. For b = 0.4 and b = 0.5, 
we find that « 90% of the dislocations in the system are 
confined to the grain boundaries whereas ~ 10% of the 
dislocations are free within the domains. We also find 
that « 5% of the dislocations within the grain boundaries 
unbind into dischnations, which occurs generally at the 
intersection of the grain boundaries. Though the free 
dischnations are absent in the system, it does not lead 
to a long range hexatic order in the system. For b = 
0.6, we find ~ 0.14, and for other values of the field 
*6 < 0.05 (for a perfect vortex lattice, = 1)- The 
small value of '^e is caused by the random orientation of 
the domains which destroys the long range orientational 
order. 

We call the intermediate field range in which the sys- 
tem breaks into regions of ordered lattice the domain 
regime (DR) . The domain regime is configurationally dis- 
tinct from the conventional picture of a disordered state 
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FIG. 2: (Color Online) The real space configuration of vor- 
tices in a region of the simulation box for A =0.01 (a), and 
0.02 (b). The = 6400. In (b), only the defects are shown. 

for which the distribution of topological defects is homo- 
geneous. In the DR, the system is quasi-ordered on the 
length scale of the domain size Rd- The vortex lattice 
shows translational and orientational order within the 
domains, even though the long range order is absent in 
the system. Fig.l(b)-(d) shows real space configurations 
in the DR as h is increased. 

For h > 0.6, the defect density increases rapidly with 
the concomitant decrease in the domain size. Small do- 
mains of ordered lattice of width 3-4ao can be seen until 
b « 0.8, as evident from Fig. 1(e). Increasing b > 0.8, 
the system becomes amorphous with an average distance 
^ flo between the defects. The defect distribution is ho- 
mogeneous throughout the system, and the configuration 
is similar to a frozen liquid. A typical real space configu- 
ration is shown in Fig. 1(f). The fraction of free disclina- 
tions is significantly higher than that observed in the DR 
but it is difficult to isolate them from the dense network 
of defects. The vortices with coordination number 4 and 
8 accounts for «6-8% of the total defects. Overall, the 
real space images in Fig.l suggests a reentrant change in 
the configuration, from a low-field amorphous to an in- 
termediate field quasi-ordered state, which finally crosses 
over to a high-field amorphous state. 

In the DR, the average domain size Rd is dependent on 
b and A. Rd ^ 5 — Tuq for low fields and increases in the 
intermediate field range. For b ~ 0.6 and A — 0.02, the 
size of some domains exceeds 20ao. By decreasing the 
pinning strength A to 0.01, we find a remarkably well 
ordered lattice with no topological defects for the sys- 
tem size = 6400, as shown in Fig. 2(a). This suggests 
that for sufhciently weak pinning strength, large domains 
of ordered lattice, comparable to sample size in typical 
experiments, can exist in 2D. Fig. 2(a) should be com- 
pared with Fig. 2(b) which shows the meandering grain 
boundaries formed by the defects for A = 0.02. With 
increasing A, Rd decreases from « 20ao for A = 0.02 
to « 3-5ao for A = 0.075. This is shown in Fig.3. For 
strong pinning, the dislocations tend to cluster in some 
regions implying that individual pinning centers locally 
tear the vortex lattice. With increasing A, the field range 
over which the DR exists is reduced and in the strong 




FIG. 3: (Color Online) The vortex configuration in a region of 
the simulation box for b — 0.6 and for A = 0.03 (a), 0.04 (b), 
0.05 (c), and 0.075 (d). The number of vortices = 4096. 



pinning limit, the system is amorphous for all values of 
b. 

In recent experiments, the changes in the real space 
configuration of vortices were studied in weak pinning 
NbSe2 samples across the order-disorder transition by 
the Bitter decoration technique jj^L l2^ . The order- 
disorder transition was previously identified in transport 
measurements and have been speculated to underlie the 
peak effect in the critical current density '27']. The deco- 
ration images show that the vortices form large ordered 
domains. The domains are separated by domain walls, 
which are defined by chains of dislocations. This domain 
formation is present throughout the B-T plane (below 
the melting line), hence the authors summarized their 
findings as the "absence of amorphous vortex matter" . 
Fasano et al. found that «85-90% of the defects are in 
the grain boundaries, whereas the remaining defects are 
isolated dislocations. All of these findings are consistent 
with our numerical findings and estimates in the inter- 
mediate field range for A = 0.02. 



Domain size distribution 

A useful quantity to characterize the DR is the distri- 
bution of the domain size N{sd), where the area of the 
domain Sd is in units of ag. Unlike in lattice models, ex- 
tracting N{sd) in models with continuous symmetry is 
not straightforward. The lattice vectors can change con- 
tinuously from domain to neighboring domain without 
nucleating defects, which makes it difficult to define the 
domain wall. In many cases, the domain walls, which are 
formed by the grain boundaries are not closed. Analy- 
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FIG. 4: (Color Online) The domain walls in the real space 
configuration, (a) and (b) shows small angle grain bound- 
aries, whereas large angle grain boundaries can be seen in (c) 
and (d). The magnetic field b = 0.50 for (a) and (d), and 
0.65 for (b) and (c). 
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FIG. 5: The histogram plot of the domain size distribution 
N(sd), where Sd is in the units of Oq. The values of b are 
(a) 0.3, (b) 0.5, (c) 0.7, and (d) 0.8. The = 4096 and 
A = 0.02. 



sis of the real space configuration suggests that the do- 
main walls are generally composed of two types of grain 
boundaries, depending on the misorientation angle Od 
between the neighboring domains. For the small angle 
grain boundaries, dd ^ 10°-16°, and the dislocations are 
separated by 3-5ao. In large angle grain boundaries, the 
dislocations form closely packed string-like structures, 
and 9d > 20° . Typical domains and domain walls formed 
by the grain boundaries are shown in Fig. 4. 

To extract N{sd), we used the following procedure to 
define domain walls in regions where the dislocations are 
apart. The mis-orientation angle 9d is obtained between 
successive vortices along one of the lattice vectors con- 
necting neighboring domains. If 9d is between 12° and 
18°, that vortex is considered as part of the domain wall. 
With this method, the domain boundary in many cases 
could be determined with reasonable accuracy. This pro- 
cedure is intended to be instructive rather than decisive, 
as it contains some arbitrariness. For example, in some 
regions the domain walls are wider than oq and then 
the mis-orientation angle is split across the domain wall. 
Also, this method does not count the really small angle 
domain walls, the ones with 9 < 12°. 

The area of the enclosed domains is used in creating 
the N{sd) histogram. Fig. 5 shows the histogram plot of 
N{sd) for various values of the magnetic field b. The 
total number of vortices Ny = 4096 and A = 0.02. At 
small and large fields, the histogram can be adequately 
characterized by a single parameter, e.g. its half-width. 
The distribution is relatively narrow with few large do- 
mains. For the intermediate fields, the N{sd) exhibits 
a broad distribution with substantial weight toward the 
tail region. This suggests that more than one parameter 
is required to characterize these distributions, especially 



the excess weight in the tail region. 



Roughness exponent 

The interaction of the vortex lattice with the quenched 
impurities leads to displacement of the vortices from 
their perfect lattice position. An important quantity 
which characterize the change in the position of the vor- 
tices is the relative displacement correlation, defined as 



iy(r) = [u(r)-u(0)f, (2) 

where the overbar represents the average over quenched 
impurities. The u(r) is the displacement of the vortex 
relative to its position in the perfect lattice. The po- 
sitional order parameter correlations Cdr) can be ex- 
pressed in terms of W{r) as Cair) ~ e~'^ w{r)/2^ where 
G is one of the reciprocal lattice vectors!^. For the 
crystalline state, W{r) = and CG{r) = 1. The ef- 
fect of the quenched impurities is to increase W{r) and 
hence reduce the positional order parameter correlations 
of the lattice. The structure factor at G, measured in the 
neutron scattering experiments, is related to the Fourier 
transform of Cq (''') . 

The roughness of an elastic medium is parameterized 
by the exponent ^, which is defined as W{r) ^ r^^ . In 
the flat phase of the medium C < 0, and in the rough 
phase C > (the C = gives logarithmic roughening with 
W{r) ~ Inr). For a 2D vortex system, there are three 
length scales which emerge in various theories depending 
upon the displacement u{r): 
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(1) r < Re- In the collective pinning theory jgj Rc 
represents the size of the region in which the vortex lat- 
tice is coherently pinned by the impurities. More pre- 
cisely, Rc is the length scale at which the displacement 
u{r = Rc) ~ ^. Rc is obtained by minimizing the total 
energy (elastic energy -I- pinning energy) and is given by 



Rc 



r 1/2 ■ 



(3) 



The 6*66 is the shear modulus of the vortex lattice and 
the average pinning force / ^ ^/fpi where rp is the 
range of the pinning potential as defined in Section II. 
For the Ki^{r/\) potential, the field dependence of the 
shear modulus h;ve been derivedH and IS given as 
Cee ~ {t^\)'i (1 ~ ^Y- III dimensionless units, the Rc 
becomes 



Rc 



,63/2(1 _ 6)2 



(4) 



The Rc increases with b and attains a maximum before 
decreasing as 6 ^ 1. 

The r < Rc regime is often referred as the random 
force (RF) regime. The roughness exponent in this 
regime is given by C = ^-5^ for a Z?-dimensional sys- 
tem. Thus, CflF = 1 and Cr^ = 0-5 for 2D and 3D 
systems, respectively. 

(2) Rc < r < Ra- Beyond Rc, the displacement u{r) 
continues to grow but with smaller exponent. Ra de- 
fines the length scale at which the positional correlation 
begins to decay, i.e., the displacement u{r — Ra) 
qq. Between Rc and Ra, the system is in the ran- 
dom manifold (RM) regime. In this regime, the rough- 
ness exponent have been obtained using a Flory type 
argument [33, HI] which gives CVai — h- A more re- 
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fined scaling argument |39j gives — 0.4. For weak 
pinning, the length scales Ra and Rc are related by 
Ra^Rciao/^)^/'^'"^'. 

(3) Ra < r < ^D- Beyond Ra, W{r) grows as|l 
W{r) ~ In^(r) as derived through a variational approach 
and confirmed by replica symmetric R,G|4lj. assuming 
the lack of dislocations at these scales [ioj. This growth 
form holds up to the length scale ^d, at which unbound 
dislocations appear. For weak pinning, ^ J?^[l8l|: 



£,D ^ Ra exp 




(5) 



where c is a temperature dependent numerical constant 
and (Tg is the impurity strength. For Ra ^ ag and low 
temperatures, can become exceedingly large and the 
system appears similar to the BG phase in 3D. 

(4) < r: Beyond ^_d, unbound dislocations lead to 
exponential decay of the positional correlation and the 
system is disordered. 
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FIG. 6: The W(r)/ag as a function of r/ao for A = 0.01 for 
iV„ = 6400. The inset shows the same for A^^ = 4096. The 
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where the exponent is shown for the random 



force (RF) regime and the random manifold (RM) regime. 



We have obtained the length scale and the roughness 
exponent of the vortex lattice in the DR. The relative 
displacement correlation W{r) was calculated using the 
following procedure. First, a crystalline state with the 
lattice constant corresponding to the value of b is con- 
structed using one of the vortex coordinates rg in the 
real lattice as the origin. The mean square displacement 
between the perfect lattice and the underlying real lat- 
tice is then minimized by varying the orientation of the 
perfect lattice relative to the real lattice. The u{r) is 
then computed for each of the vortices. This procedure 
is repeated for different rg's, and the W{r) is computed 
by averaging over all rg's. 

Fig. 6 shows the plot of VF(r)/ag for b — 0.6 and A = 
0.01 for iV„ = 6400. W{r) for iV„ = 4096 is shown in 
the inset. For these parameters, even the largest Ny = 
6400 system is free of dislocations (see Fig. 2(a)). Due to 
periodic boundary condition, the length scale probed in 
the simulation is half the system size, which for = 
6400 and 4096 are 40ag and 32ag, respectively. W{r) 
shows an initial power-law increase with an exponent 
( w 0.65 — 0.72. The system exhibits a crossover around 
r* « 4 — 5ag into a regime where the increase of W{r) 
is slower. Between 5ag ^ r < 15ag the exponent is 
C « 0.40 - 0.42. For r > 16ag, the growth of W{r) slows 
down considerably. 

It is plausible that r* = Rc, where the system crosses 
over from the RF regime to the RM regime. This can 
be verified by calculating Rc using Eq.(4), which gives 
a value « 4ag for b = 0.6, in reasonable agreement with 
the value of r*. In our system ^ = O.lAg, and at & = 0.6 
the lattice constant ag = 0.347Ag. From Fig. 6, u{Rc) « 
0.2ag « 0.07Ag, thus confirming that u{r = Rc) ~ ^■ 

Beyond Rc, the length scale Ra of the RM regime 
is defined as u(r — Ra) ^ ag. This translates to 
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W{r) « 0.12ag for b = 0.6. From Fig.6, we find that 
wlr) flattens at O.f at r - 13ao for = 6400, which 
suggests that Ra « 13ao (for Ny — 4096, Ra ~ fSao). 
Beyond Ra, the growth of W{r) slows down consider- 
ably which indicates the appearance of the asymptotic 
regime. Within the qBG theory W(r) is expected 
to grow as In^(r) in the asymptotic regime. This be- 
havior unfortunately could not be verified due to insuf- 
ficient range of data points. In sum, we identify the 
r ~ IflQ — 5ao as the RF regime, the r ^ 5ao — 15ao as 
the RM regime, and r > 15ao as the asymptotic regime. 

The value of C ~ 0.65 — 0.72 obtained from the sim- 
ulation in the RF regime is smaller than the theoretical 
prediction for — 1- We speculate that the interac- 
tion Ko{r/X) between the vortices in 2D increases the 
stiffening of the vortex lattice at short distances which 
leads to weaker roughening. In the RM regime, the ex- 
ponent C « 0.40 — 0.42 is in good agreement with the 
value of 0.4 expected from the scaling argumentjUj. Us- 
ing = 0.4, the value of Ra « 15oo is much smaller 
than the value Ra ~ RciaQ/S^Y^^R^" . This is possibly 
related to the large magnetic field h — 0.6 used in ob- 
taining W{r). For this field, ao is comparable to ^, and 
Rc is large compared to smaller magnetic fields. The 
RM regime is expected to disappear for oq ~ S,, and 
have been shown in the case of 3D system |42|. 

An interesting outcome of the above analysis of W{r) 
is that the average domain size Rd ^ Rc, and hence, the 
collective pinning theory cannot account for the appear- 
ance of domains in the intermediate fields. The asymp- 
totic regime in W{r) suggests that the qBG theory is 
qualitatively correct. Within the qBG theory, the dis- 
tribution of dislocations beyond the length scale is 
expected to be homogeneous, unlike the grain bound- 
ary formation observed in our simulation. One possi- 
ble way to account for the grain boundary formation 
is to consider the long range interaction between the 
dislocations. Since the interaction between the disloca- 
tions is anisotropic, for some values of dislocation den- 
sity, the grain boundaries may lead to a lower energy 
state. This is also supported from a recent work on 
dipole svstems|4^. At low densities, the dipoles exhibit a 
gaseous phase, and their distribution is roughly homoge- 
neous. At higher densities the phase is characterized by 
dipoles forming chains or strings. Since the dislocations 
of the vortex lattice are in fact dipoles of disclinations, 
these results are quite analogous to our identification of 
a domain regime in the vortex matter. 

Defeat density 

The three field regimes discussed in the context of the 
real space configuration can also be inferred from the be- 
havior of the defect density nd{h). Fig. 7(a) shows nd{h) 
for A = 0.02 and = 4096. The behavior for smaller 




FIG. 7: (a) The topological defect density nd{h) for A = 0.02 
and A'„ = 4096. Also shown is the nd{b) for the smaller 
system size, (b) The fdi^) for b = 0.6. Inset: The plot of fd 
as a function of A'^^ for A = 0.02 and b — 0.6. There is a 
critical system size below which dislocations are not present 
in the system. 



system sizes (Af^=800-1200) is also shown on the same 
plot. The nd(b) increases linearly in the low field amor- 
phous regime. Above a crossover field 6/ « 0.1, nd{b) 
flattens and becomes weakly field dependent in the DR. 
For 6 > 0.6, nd{b) increases rapidly, and above bh ~ 0.8 
the system crosses over to the high field amorphous 
regime. It is possible to define a length scale Ld ~ i^d^^"^, 
as the nominal average defect separation. For 6 = 0.6 
(domain regime) Ld ~ Sag, which is much smaller than 
even Rc and does not correspond to any feature in the 
real space configuration, and reflects the highly inhomo- 
geneous nature of the defect distribution in the domain 
regime. On the other hand, in the high-field amorphous 
regime Ld ^ ao, which is also the distance between the 
defects, thus reflecting homogeneity of the distribution 
of defects. 

The nd{b) in Fig. 7(a) shows strong similarity with 
the experimental observation in 2D system of magnetic 
bubbles 44J. In Ref. ,,44] . the intermediate regime was in- 
terpreted as the hexatic phase and the high field amor- 
phous regime as the isotropic hquid phase^^. Later 
simulation^^ also suggested a T = dislocation unbind- 
ing transition driven by disorder. As discussed above, 
the presence of domain walls (grain boundaries) in our 
system suppresses the long range orientational order. 
This rules out the possibility of a transition between 
the hexatic phase and the isotropic liquid phase as the 
underlying reason for the rapid increase in nd{b). How- 
ever, a rapid crossover, similar to that predicted between 
the qBG at low temperatures and vortex liquid at high 
temperaturesfigj . is still possible between the DR and 
the high field amorphous regime, especially at weaker 
pinning where the domain size Rd is large|47j. 

For smaller system size {Ny = 800 — 1200), a topo- 
logically ordered phase appears in the intermediate field 
range in which = (see Fig. 7(a)). This is a finite 
size effect, which reflects the sensitivity of the DR to 
the system size L. For L < (^d, the DR can appear as 
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FIG. 8: (Color Online) The domains in a region of the sim- 
ulation box for (a) b = 0.7 and (b) 0.8. The = 4096 and 
A = 0.02. 

a topologically ordered state free of dislocations. This 
implies that for a given system size, there exists a criti- 
cal Ac below which dislocations are not favored. This is 
observed in the simulation, as shown in Fig. 7(b) where 
the defect fraction /d(A) (number of defects per vortex) 
is plotted for b = 0.6. For the smaller system Ny = 900, 
fd goes to zero at Ac — 0.03, and the system exhibits 
an ordered phase for A < Ac. Increasing iV„ to 4096 re- 
duces Ac to 0.01, and in the asymptotic limit Ny — > oo 
(hence, L oo), we expect Ac 0. In the DR, fd does 

1/2 

not increase continuously with increasing i(cx K'") but 
shows a sharp jump from the dislocation-free state to the 
domain state at a characteristic system size as shown in 
the inset of Fig. 7(b). 

Crossover from DR to high-field amorphous regime 

As discussed above, the nd{b) shows a rapid crossover 
from the DR to the high field amorphous regime. To 
understand the mechanism for this sharp crossover, we 
identified the domains and the domain walls between 
b = 0.5 and b = 0.8 for A = 0.02. For the intermedi- 
ate fields b « 0.5-0.6, the grain boundaries are generally 
smooth and nd{b) is weakly field dependent. For b > 0.6, 
the rapid increase in nd{b) occurs within the domain 
walls. Consequently, the domain wall length increases, 
which is accommodated through enhanced roughening of 
the domain walls. This is evident from Fig.8(b). The in- 
crease in the roughening also facilitates the unbinding of 
the dislocations into free disclinations and subsequently 
drives the crossover into the VG state. In such a sce- 
nario, we conjecture that domain walls undergo disorder 
driven roughening transition at the crossover between 
the DR and the high field VG. It would be of interest to 
obtain the domain wall roughening exponent across the 
crossover regime. 

Finite temperature simulation 

In this section, we compare the current anneal- 
ing method with the conventional simulated annealing 
method, as it is well known that different sample prepa- 




FIG. 9: (Color Online) (a) The defect fraction fd{T) for var- 
ious magnetic fields b with A = 0.03. The temperature T 
is decreased slowly from the high-T liquid phase, (b) The 
fd{b) at T = obtained from thermal annealing and current 
annealing. 

ration techniques can result in the vortex system not 
reaching its equilibrium configuration '33", '48] . In the lat- 
ter method, the temperature T is reduced from the high 
temperature liquid phase slowly so as reach thermal equi- 
librium at each value of T. This method is commonly 
used to search for the ground state of disordered systems. 

For the thermal annealing, the system was equili- 
brated for 5 X 10^ — 1 X 10^ time steps before averaging 
over a similar time scale to calculate the defect fraction 
fdiT). The number of vortices N., = 900 and A = 0.03. 
For this system size, the CA method gives a topologi- 
cally ordered phase for b between 0.6 and 0.75. Fig. 9(a) 
shows fd{T) for various values of the magnetic field. As 
the temperature is lowered, for b — 0.2 fd{T) decreases 
monotonically to a finite value with ^ slowly varying. 
There is no evidence of a transition as a function of the 
temperature. With increasing 6, the slow freezing is re- 
placed by a sharp decrease in fdiT) at a particular tem- 
perature Tm, similar to the equilibrium melting transi- 
tion. For b = 0.65, fdiT) at decreases by « 76% 
of the value above Tm- For b > 0.8, the melting-like 
transition is again replaced by slow freezing of the high 
temperature liquid phase. 

In Fig. 9(b), fdib) at T = obtained by TA is com- 
pared with that obtained by CA. At intermediate fields, 
the thermally annealed samples exhibit the presence of 
dislocations already at these smaller systems sizes. As 
described above, the current annealing method requires 
larger systems sizes to correctly display this same phe- 
nomenon. Otherwise, the two curves track each other 
very closely over most of the field range, including the 
low field slow decay of fdib) and the rapid rise at high 
fields. For the intermediate fields, the T = configura- 
tion obtained from TA also shows grain boundary for- 
mation, similar to that observed from the CA method. 



CONCLUSION 

We have presented a detailed numerical analysis of the 
real space configuration of 2D vortex system in the pres- 
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ence of quenched impurities. For weak pinning, the dis- 
ordered state in the intermediate field range is inhomoge- 
neous. The majority of the dislocations in this state are 
confined to grain boundaries, which form domain walls 
between regions of topologically ordered vortex lattice. 
There are no free disclinations in the system. This state 
is referred as the domain state and the intermediate field 
range as the domain regime. 

The domain size distribution N[sd) was calculated in 
the domain regime. N{sd) shows a broad distribution 
with a large weight in the tail region at intermediate 
fields. Therefore, more than one length scale is required 
to properly characterize the domain size distribution in 
the domain regime. With increasing b, the distribution 
becomes narrow and the peak shifts toward the origin. 
For weak pinning, the size of the domains can become 
exceedingly large. 

The domain regime is bounded by an amorphous 
regime at low fields and high fields. The defects in the 
amorphous regime are separated by the smallest length 
scale ~ oo and show homogeneous distribution unlike the 
grain boundary formation in the domain regime. The 
domain regime shows rapid crossover into the high field 
amorphous regime. From the changes in the configura- 
tion, we identified the roughening of the domain walls as 
the plausible mechanism driving the rapid crossover. 

The relative displacement correlation W{r) in the do- 
main state was also calculated for weak pinning. Three 
distinct regimes were observed: a random force regime, 
a random manifold regime and the asymptotic regime. 
Crossover from random force regime to the random man- 
ifold regime is found to occur at Rc ^ 4 — Soq. The value 
of Rc agrees with that obtained from the collective pin- 
ning theory. The roughness exponent in the random 
manifold regime is found to « 0.40, within the range of 
various theoretical predictions. 

The observation of random manifold and asymptotic 
regimes within the domains for weak pinning suggests 
that the vortex lattice is correctly described by the qBG 
idea, though the exact form of the W{r) could not be 
ascertained. At length scales greater than the domain 
size, the appearance of the domain wall formed by dis- 
locations is not captured by the quasi-Bragg glass the- 
ory. Therefore, it remains to be seen whether besides the 
domain regime the 2D vortex matter supports a quasi- 
Bragg glass where the dislocations are homogeneously 
distributed. 
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